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I.  Introduction 


The  representation  of  antenna  and  scattering  patterns  in  three  dimen¬ 
sions  provides  a  useful  tool  for  analyzing  power  flow  or  field  strength,  es¬ 
pecially  when  used  in  conjunction  with  field  line  contour  plots  (Ref.  1). 
Geometrically,  these  patterns  represent  surfaces.  A  comprehensive  bibliog¬ 
raphy  of  surface-representation  algorithms  is  given  by  Griffiths  (Ref.  2). 
Most  of  these  algorithms  are  geared  towards  representing  the  complex  shape 
of  a  physical  object.  Specific  algorithms  for  plotting  mathematical  functions 
of  two  variables  in  rectangular  coordinates  generally  follow  the  method  pro¬ 
posed  by  Wright  (Ref.  3).  Advanced  line  drawing  algorithms  suitable  for 
cylindrical  and  spherical  coordinates  have  been  developed  by  Scott  (Ref.  4). 
Unfortunately,  all  these  methods  require  a  significant  amount  of  additional 
computation  to  be  able  to  plot  the  surface  with  the  hidden  lines  removed. 
The  algorithm  developed  and  presented  in  this  paper  avoids  this  extra  com¬ 
putation  by  exploiting  known  properties  of  the  surfaces  being  plotted. 

In  general,  line  drawing  algorithms  do  not  remove  hidden  lines,  but 
rather,  simply  do  not  draw  them.  An  alternative  way  to  remove  the  hid¬ 
den  lines  of  a  surface  is  to  paint  over  the  hidden  part  with  the  same  color  as 
the  background.  This  is  exactly  the  technique  an  artist  would  use  to  paint 
a  landscape.  The  image  is  placed  on  the  viewing  surface  from  background 
to  foreground  and  hidden  lines  are  painted  over.  Of  course,  this  requires  the 
graphics  system  to  be  able  to  fill  or  erase  a  polygonal  region.  Therefore, 
the  algorithm  proposed  in  this  paper  is  not  suitable  for  representing  surfaces 
by  means  of  mechanical  pen  plotters.  However,  it  is  ideally  suited  for  video 
displays  and  laser  printers. 

The  algorithm  is  based  on  the  following  postulate.  If  a  function  /(«,  t>), 
where  u  and  v  are  two  coordinates  of  an  orthogonal  system,  generates  a 
single- valued  surface  in  the  variables  u  and  v,  then,  there  exists  a  systematic, 
although  not  unique,  ordered  sequence  in  which  to  draw  the  surface  from 
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back  to  front.  This  sequence  is  known,  a  priori,  once  tho  observation  angles 
are  specified.  Therefore,  no  hidden  line  removal  is  necessary.  Thus,  plotting 
a  function  with  the  hidden  lines  removed  takes  the  same  amount  of  time  as 
plotting  without  removing  the  hidden  lines.  A  single- valued  surface  is  defined 
as  a  surface  in  which  there  is  a  one  to  one  correspondence  between  each 
pair  of  coordinates  (u,  u)  and  a  point  on  the  surface.  Aperture  distributions, 
antenna  patterns,  and  scattering  patterns  are  known  to  generate  surfaces 
that  are  singi°- valued  with  respect  to  a  particular  coordinate  system. 


II.  Coordinate  Systems 


The  functions  of  interest  in  this  study  are  plotted  in  rectangular,  cylin¬ 
drical,  and  spherical  coordinates.  The  standard  variables  that  describe  these 
coordinate  systems  are  defined  by 


r  —  yjx2  +y2 

(1) 

R  =  yjr2  +  z2 

(2) 

<f>  =  tan-1  (y/x) 

(3) 

6  =  tan-1  (r/z)  - 

(4) 

In  order  to  organize  the  algorithm  inputs  consistently,  let  f(u,  v)  describe  a 
function  of  two  variables  from  one  of  these  coordinate  systems.  Table  1  shows 
the  convention  adopted  for  the  relationship  between  u  and  v  and  the  stan¬ 
dard  variables.  Since  the  final  plot  is  placed  on  a  two-dimensional  surface, 
it  is  convenient  to  perform  this  step  in  rectangular  coordinates.  Therefore, 
all  plotting  points  will  be  converted  to  rectangular  coordinates  prior  to  any 
drawing.  The  conversion  conventions  are  shown  in  Table  2.  Once  the  points 
on  t>e  surface  are  converted  to  rectangular  coordinates,  the  graphical  oper¬ 
ations  of  scale,  rotation,  and  projection  may  be  applied. 

Table  1.  Variable  Definitions 
cylindrical 

u  as  <j>  <f> 

v  y  t  z 

Table  2.  Conversions  to  Rectangular  Coordinates 


rectangular 


spherical 


rectangular 


cylindrical 


spherical 


III.  Scale,  Rotation,  and  Projection 

The  surface  defined  by  f(u,v)  is  represented  by  the  graphical  transfor¬ 
mations  of  scale,  rotation  and  projection.  Each  graphical  operation  must  be 
applied  to  individual  coordinates  of  the.  surface.  Let  each  point  of  the  surface 
be  defined  by  the  vector  w  such  that 


The  graphical  operations  may  now  be  defined  as  matrices  that  operate  on 
this  vector. 

In  order  to  enhance  some  visual  attributes,  it  may  be  desirable  to  scale 
each  point  before  plotting.  A  scaling  matrix  S  is  defined  as 


where  each  component  is  some  specified  constant.  Note  that  to  preserve  the 
linearity  of  the  scale,  the  condition  S«  =  Sy  must  be  satisfied  in  cylindrical 
coordinates  and  Sx  =  Sy  =  Sz  must  be  satisfied  in  spherical  coordinates.  The 
coordinate  points  are  then  scaled  by  forming  the  matrix  vector  product  Sw. 
Another  type  of  scaling,  which  is  quite  commonly  used  in  plotting  antenna 
and  scattering  patterns,  involves  converting  the  function  to  decibels.  This 
allows  the  viewer  to  observe  more  detail  of  the  sidelobe  bel  avior.  Let  /  be 
normalized  to  the  range  0  <  /  <  1.  Then  define  the  zero  reference  in  dB  as  v. 
Next,  define  a  plot  floor  level  in  dB  as  tj,  such  that  77  >  v.  This  plot  floor  is 
the  level  to  which  the  function  is  set  to  for  any  value  below  the  floor.  This 
avoids  a  cluttered  graph  that  results  from  too  many  low-level  sidelobes.  The 
function  /  can  then  be  converted  to  a  dB  scale  such  that  0  <  /dB  <  1  by  the 
nonlinear  transformation 


/dB  = 


/<lOW10) 
j£y(10log(/)  +  M)  />  10<”/10). 


(7) 
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The  observer  is  assumed  to  be  stationary,  so  the  surface  of  the  function 
must  be  rotated  to  the  correct  view.  Viewing  any  finite  three-dimensional 
object  requires  a  minimum  of  one  rotation  axis.  However,  it  is  usually  neces¬ 
sary  to  have  rotation  around  two  axes.  These  axes  should  be  perpendicular 
to  allow  the  widest  range  of  viewing  angles.  The  convention  adopted  in  this 
study  is  to  allowr  an  azimuthal  rotation  around  the  z  axis  and  an  elevation 
rotation  in  the  yz  plane.  This  is  conveninent  for  plotting  antenna  patterns, 
aperture  distributions  and  scattering  patterns.  The  visual  result  is  a  graph 
that  appears  to  spin  in  azimuth  around  the  z  axis  and  is  tipped  in  elevation 
toward  or  away  from  the  viewer. 

It  is  important  to  distinguish  between  the  observation  angles  and  the 
rotation  angles.  The  viewer-supplied  elevation  observation  angle  is  defined 
as  0O j  and  the  azimuth  observation  angle  is  defined  as  <£0.  These  angles  will 
be  restricted  to  the  ranges  0  <  0O  <  *  and  0  <  <j>0  <  2*-,  respectively.  The  use  of 
these  angles  provides  viewers  a  familiar  frame  of  reference.  In  order  to  view 
the  surface  from  these  angles,  it  is  necessary  to  define  an  azimuth  rotation 
angle  a  and  an  elevation  rotation  angle  /?.  These  angles  are  dependent  on 
the  observation  angles  and  the  orientation  of  the  viewer  to  the  rectangular 
coordinate  system.  It  is  assumed  that  the  observer  will  look  in  from  the 
positive  z  axis  onto  the  zy  plane.  The  *  axis  increases  from  left  to  right  and 
the  y  axis  increases  from  bottom  to  top.  Once  this  convention  is  established 
and  the  observation  angles  are  specified,  the  rotation  angles  can  be  calculated 
as 


a~"~4>o  (8) 

/?  =  90.  (9) 

The  order  of  rotation  is  not  commutative.  The  plots  must  first  be  spun  in 
azimuth  by  the  angle  a  and  then  tipped  in  elevation  by  the  angle  0.  Reversing 
the  order  would  allow  the  plot  to  tip  from  side  to  side,  which  presents  an 
awkward  picture.  The  azimuth  rotation  in  the  zy  plane  is  represented  by  the 
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(10) 


matrix  vector  product  RQw,  where 


R  a 


cos  a  —  sin  a  0 
sina  cosa  0 
0  0  1 


and  the  elevation  rotation  in  the  yz  plane  by  the  product  R^w,  where 


R fi  = 


l 

0 

0 


0  0 
cos/3  sin/3 

—  sin/3  cos/3 


(11) 


The  complete  transformation  can  be  represented  by  the  matrix  T,  defined 
by 


T  =  R^RoS.  (12) 

If  w'  represents  the  transformed  coordinates  then 

w'  =  Tw.  (13) 

Viewing  a  three-dimensional  plot  on  a  two-dimensional  surface  requires  a 
projection  of  the  three-dimensional  coordinate  points  onto  a  two-dimensional 
viewing  plane.  Perspective  projection  implies  that  the  further  a  point  is 
away  from  the  viewer,  the  smaller  it  appears.  This  type  of  projection  is 
appropriate  for  images  that  evoke  a  strong  depth  cue  such  as  a  building 
or  a  landscape.  However,  mathematical  functions  do  not  require  a  strong 
depth  cue  since  there  is  no  physical  object  being  represented.  Therefore,  it 
is  sufficient  to  project  each  point  along  a  parallel  line  until  it  intercepts  the 
viewing  plane.  This  type  of  projection  is  called  parallel  projection.  Using 
parallel  projection  implies  that  only  the  s'  and  j/  components  of  the  vector  w 
are  needed  to  represent  the  surface.  The  z'  coordinate  represents  the  depth 
of  each  point  and  is  not  used. 

These  three  graphical  transformations  are  independent  of  the  type  of 
function  being  plotted.  They  represent  the  transformation  of  a  point  on 
the  function  surface  to  a  point  on  the  viewing  surface.  The  order  in  which 
points  arc  operated  on  by  the  transformation  is  determined  by  the  sequencing 
algorithm  discussed  in  the  nex*  section. 
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IV.  Algorithm  Description 


The  proposed  algorithm  is  based  on  the  generic  “painter’s  algorithm.” 
This  concept  means  that  the  surface  of  the  plot  is  rendered  by  building  up  the 
image  from  back  to  front.  Since  the  function  has  two  variables,  the  surface  is 
most  naturally  described  by  a  collection  of  quadrilaterals.  Each  quadrilateral 
is  placed  on  the  viewing  surface  and  is  painted  with  the  background  color  of 
the  surface.  A  line  around  the  perimeter  is  then  drawn.  The  key,  of  course, 
is  to  know  the  order  in  which  to  place  the  quadrilaterals  on  the  viewing 
surface.  A  previous  method,  developed  by  the  author  (Ref.  5),  used  the 
transformed  z  coordinate  of  the  centroid  of  each  quadrilateral  and  sorted  by 
depth.  This  method  was  general  enough  to  plot  any  surface.  However,  the 
time  required  to  compute  the  depth  of  each  centroid  and  to  sort  the  values 
led  to  a  significant  time  delay.  The  new  method  proposed  avoids  this  delay 
by  drawing  in  a  prescribed  order  based  on  the  observation  angles  <f>0  and  60. 

The  variables  u  and  v  are  discretized  such  that  u  =  um  for  m  =  1  ...JV„  and 
v  =  t>„  for  n  =  1  ...Nv.  By  convention,  the  values  sequence  from  the  minimum 
to  the  maximum  values.  The  function  /( u,v)  is  then  sampled  at  the  discrete 
points  fm,n  =  /(%,«! i).  Each  quadrilateral  has  a  reference  corner  that  has 
index  values  (m,  n).  The  other  three  corners  are  dependent  on  these  indices 
and  are  given  by  (m-f-l.n),  (m+l,»i+l),  and  (m,n+l).  It  should  be  noted  that 
the  definition  of  a  quadrilateral  is  extended  to  allow  any  number  of  corners 
to  have  the  same  location.  Therefore,  a  point,  a  line,  and  a  triangle,  can 
also  be  represented  by  a  quadrilateral.  The  crux  of  the  algorithm  is  to  find 
a  systematic  method  of  sequencing  through  the  indices  in  order  to  approxi¬ 
mate  a  back-to-front  ordering.  The  sequencing  orders  presented  for  the  three 
coordinate  systems  analyzed  were  chosen  based  on  ease  of  programming. 
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V.  Rectangular  Coordinates 


Rectangular  coordinates  are  the  most  common  way  of  representing  pat¬ 
terns  in  three  dimensions.  The  coordinates  *  and  y  are  replaced  by  <p  and 
0.  Assume  for  simplicity  that  the  coordinates  are  translated  so  that  the 
origin  is  in  the  interior  of  the  range.  This  does  not  affect  the  generality  of 
the  algorithm,  but  merely  the  presentation.  Since  the  function  is  plotted 
on  a  rectangular  base,  it  may  be  surmised  that  one  corner  of  the  surface 
will  always  be  nearest  to  the  observer  and  the  opposite  corner  will  be  the 
farthest  away.  Because  of  the  rotation  conventions  used,  this  is  dependent 
only  on  the  azimuth  observation  angle  4>q  and  independent  of  the  elevation 
observation  angle  90.  Therefore,  the  first  part  of  the  algorithm  determines 
which  corner  is  nearest  to  the  observer.  Once  the  orientation  is  established, 
the  quadrilaterals  are  simply  drawn  from  the  back  corner  to  the  front  cor¬ 
ner  by  rows  that  alternate  in  direction.  Alternating  the  rows  helps  to  move 
forward  in  a  more  uniform  manner.  Figure  1  shows  the  drawing  directions 
as  the  quadrilaterals  are  placed  on  the  viewing  surface  from  back  to  front. 
Note  that  at  the  angles  0°,  90°,  160°,  270°  and  360°  there  are  two  back  corners 
equidistant  from  the  viewer.  In  this  case,  it  does  not  matter  which  corner  is 
chosen  as  long  as  it  is  chosen  consistently.  Figure  2  shows  the  far-field  power 
pattern  of  a  uniformly  excited  square  aperture  plotted  at  observation  angles 
(f>0  =  30°  and  60  =  60°.  The  plotting  algorithm  used  is  given  in  pseudo-code  as 
follows: 

M  —  minimum  of  {Nu,  Nv) 
if  0  <  <j)0  <  tt/2  then  {find  nearest  corner} 
i  =  1,  j  =  l 

else  if  7r/2  <  <f>o  <  ft  then 

*  =  -1,  3  =  1 

else  if  x  <  <fr0  <  37r/2  then 
i  =  -1,  j  =  -1 
else  if  3x/2  <  4>o  <  2k  then 

*  =  1,  3  =  -1 

end  if 

loop  from  /  =  1  to  M  -  1 

if  0  <  <po  <  7r/2  then  {find  nearest  corner) 
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Fig.  1.  Drawing  flow  pattern  for  rectangular  coordinates 
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Fig.  2.  Power  pattern  of  uniformly  excited  square  aperture  viewed  from  0O  =  60° 
and  <po  =  30° 
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mo  =  l  -  1,  n0  =  l 
mi  =  l,  nt=l 
else  if  7r/2  <  <f>o  <  ^  then 

mo  =  Nu  -  l  +  1,  no  =  l 
m\  =  Nu  -  l,  ni  =  l 
else  if  7T  <  <  37r/2  then 

mo  =  Aru  -  l  +  1,  no  =  Nv  -  l 
mi  =  Nu  -  l,  ni  =  Nv  - 1 
else  if  37t/2  <  <  2tt  then 

mo  =  l  -1,  n0  =  Nv  -  l 
mi  =  l,  ni  =  Nv  —l 
end  if 

loop  from  k  =  1  to  Nu  —  l  —  1 
m  =  mo  +  ik,  n=  no 
{  get  (um,vn)  and  other  3  corners  } 

{  convert  to  rectangular  coordinates  -  Table  2  } 
{  scale,  rotate,  and  project  using  Eqn.  (13)} 

{  fill  quadrilateral,  then  draw  perimeter  } 
continue  k  loop 
loop  from  k  =  1  to  Nv  —  l 
m—  mi,  n  —  ni  +  jk 
{  get  (wm,v„)  and  other  3  corners.} 

{  convert  to  rectangular  coordinates  -  Table  2  } 
{  scale,  rotate,  and  project  using  Eqn.  (13)} 

{  fill  quadrilateral,  then  draw  perimeter  } 

continue  k  loop 
continue  l  loop 
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VI.  Cylindrical  Coordinates 


Cylindrical  coordinate  functions  of  the  form  z  =  f(<t> ,  r)  and  r  =  f(<f>,  z)  are 
commonly  encountered  in  antenna  and  scattering  analysis.  Far-field  patterns 
can  be  plotted  in  the  coordinates  f(<f>,  r)  by  letting  r  =  6.  This  type  of  plot 
is  good  for  observing  the  finer  details  of  the  sidelobe  structure.  It  is  also  a 
natual  way  to  plot  circular  aperture  distributions.  Alternatively,  functions 
of  the  form  r  =  z)  are  useful  for  observing  near-field  patterns  or  surface 
currents  on  structures  such  as  a  cylinder  or  a  body  of  revolution. 

Cylindrical  functions  of  the  form  z  =  f{<f> ,  r)  are  plotted  on  a  circular  base 
with  a  specified  foreground  angle  4>o  and  a  background  angle  <f>a  =  <j> 0±7r,  chosen 
such  that  0  <  tpa  <  27r.  Therefore,  the  <f>  dependency  of  the  quadrilaterals 
should  be  drawn  starting  with  the  background  angle  <f>a  and  procceeding 
toward  the  foreground  observation  angle  <p0 .  The  method  employed  is  to  find 
the  index  m  of  the  um  nearest  to  (j>a  and  then  alternately  increase  and  decrease 
the  value  to  draw  from  back  to  front.  For  simplicity,  assume  that  umin  =  0 
and  umax  =  2ir.  Since  <j>  is  periodic,  the  index  m  must  be  periodic  with  period 
Nu-  1.  Drawing  the  radial  dependence  from  back  to  front  is  dependent  on 
the  value  of  <j>.  Observation  of  the  drawing  flow  pattern  in  Fig.  3  shows  that 
the  radial  dependence  should  be  drawn  from  rmax  to  rm<„  when  the  angle  um 
is  greater  than  90°  from  $0)  and  from  rm<„  to  rmax  when  the  angle  is  less  than 
90°  from  <p0-  Figure  4  shows  the  power  pattern  of  a  uniformly  excited  square 
aperture  plotted  in  dB,  with  v  =  -40  dB  and  77  =  -25  dB.  In  some  plotting 
situations  it  is  desirable  to  remove  an  angular  sector,  to  allow  better  viewing 
of  the  surface.  This  can  be  easily  accomplished  because  of  the  way  the  4> 
values  are  alternated  when  plotting.  Defining  N,  as  the  number  of  segments 
to  remove,  the  total  number  of  0  values  is  just  reduced  by  that  amount. 
Figure  5  illustrates  the  use  of  angular  cuts  in  a  plot  of  the  magnitude  of  the 
far-field  electric  field  of  a  uniformly  excited  circular  aperture.  The  following 
pseudo-code  implements  the  algorithm: 
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Fig.  3.  Drawing  flow  pattern  for  cylindrical  coordinates 
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Fig.  4.  Power  pattern  of  uniformly  excited  square  aperture  viewed  from  80  =  60° 
and  <j>o  =  30°  and  plotted  in  dB,  with  u  =  -40  dB  and  77  =  -25  dB 


A 
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Fig.  5.  Magnitude  of  far-field  electric  field  of  a  uniformly  excited  circular  aper¬ 
ture  viewed  from  Go  -  SO’  and  <&,  =  30°  and  plotted  in  dB,  with  N,  -  14, 
v  =  -40  dB  and  t)  -  -25  dB 


loop  from  k  =  1  to  Nu  -  1  {find  index  of  angle  <f>a} 
if  rpa  >  Uk  and  <f>a  <  then 

m  =  k 
end  if 

continue  k  loop 
loop  from  l  =  1  to  Nu  -  N,  -  1 
m  =  m—  (— 1)!(1  —  1) 

if  (m  <  1)  then  {force  index  to  be  periodic} 
m  =  m  +  Nu  —  1 
else  if  (m  >  Nu  -  1)  then 
m  =  m  —  (Nu  —  1) 
end  if 

if  \um  -  (j) 0|  <  ir/2  or  \um  -  cf>0 1  >  37r/2  then 

7t0  =  o,  i  =  1  {draw  from  rm{„  to  rmax} 

else 

n0  =  Nv,  i  =  -l{draw  from  rmax  to  rmilt} 
end  if 

loop  from  k  =  1  to  Nv  —  1 
n  =  no  +  ik 

{get  (um,vn)  and  other  3  corners  } 

{convert  to  rectangular  coordinates  -  Table  2  } 

{scale,  rotate,  and  project  using  Eqn.  (13)} 

{fill  quadrilateral,  then  draw  perimeter  } 
continue  k  loop 
continue  l  loop 

The  algorithm  for  plotting  functions  of  the  form  r  =  z)  uses  the  sa  me 
procedure  for  handling  the  <j>  variation.  The  z  variation  is  handled  very  easily 
by  simply  noting  that  if  the  observation  angle  60  is  less  than  or  equal  to  90°, 
then  draw  from  zmin  to  zmax.  Otherwise,  if  0o  is  greater  than  90°,  then  draw 
from  zmax  to  Zmin.  Figure  6  shows  a  section  of  circular  waveguide  and  the 
magnitude  of  the  longitudinal  surface  current  induced  on  the  inner  walls  for 
the  TM2 1  mode.  The  pseudo-code  for  the  algorithm  is  given  as  follows: 

loop  from  k  =  1  to  Nu  -  1  {find  index  of  angle  <f>a} 
if  <f>a  >  uy  and  <f>n  <  then 
ma  =  k 
end  if 

continue  k  loop 
loop  from  k  =  1  to  Nv  -  1 
if  t?o  <  7r/2  then 

n  =  k  {draw  from  zmin  to  zmax} 

else 

n  =  Nr  -  k  {draw  from  zmax  to  zmin } 
end  if 
m  =  ma 
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Fig.  6. 


Circular  waveguide  and  magnitude  of  surface  current  induced  on  inside 
wall 
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loop  from  l  =  1  to  Nu  —  N,  —  1 
m  =  m—  (—l)1  (l  —  1) 

if  ( m  <  1)  then  {force  index  to  be  periodic} 
m  =  m  +  Nu  —  1 
else  if  (m>  Nv  -  1)  then 
m  =  m  —  ( N-u  -  1) 
end  if 

{get  ( um,v„ )  and  other  3  corners  } 

{convert  to  rectangular  coordinates  -  Table  2  } 
{scale,  rotate,  and  project  using  Eqn.  (13)} 
{fill  quadrilateral,  then  draw  perimeter  } 
continue  l  loop 
continue  k  loop 
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VII.  Spherical  Coordinates 


Thu  flow  of  power  in  far-field  antenna  and  scattering  patterns  is  best 
described  by  using  plots  in  spherical  coordinates.  Establishing  the  front 
and  back  orientation  requires  knowledge  of  both  observation  angles.  The 
azimuth  background  angle  <f>a  is  defined  the  same  as  in  cylindrical  coordinates. 
Therefore,  the  <f)  variation  is  drawn  in  exactly  the  same  way  as  for  cylindrical 
plots.  The  elevation  background  angle  is  defined  as  0«  =  x  -  90.  As  shown  in 
Fig.  7,  the  elevation  dependence  is  drrwn  in  opposite  directions  away  from 
0a  when  um  is  greater  than  90°  from  <j>0.  When  Um  is  less  than  90°  from  <f>0,  the 
elevation  dependence  is  drawn  from  the  minimum  and  maximum  limits  to 
90.  Figure  8  shows  a  power  pattern  of  a  circular  aperture  plotted  in  dB  with 
v  =  t]=.  -60  dB.  As  with  cylindrical  coordinates,  it  is  sometimes  convenient  to 
cut  away  an  angular  piece  of  the  plot  so  that  the  detail  of  the  sidelobe  pattern 
may  be  observed.  However,  since  the  lobes  of  spherical  plots  are  generally 
closed  surfaces,  it  looks  better  if  the  angular  cut  is  filled  with  the  background 
color.  This  is  illustrated  in  Fig.  9.  Filling  this  cut  requires  making  a  polygon 
using  all  the  0  values  at  a  constant  <j> .  The  following  pseudo-code  implements 
the  algorithm: 

loop  from  k  =  1  to  Nu  -  1  {find  index  of  angles 
if  <f>0  >  uj,  and  4>q  <  u*+I  then 
ttiq  =  k 
end  if 

if  $<t  >  «»t  and  <f>a  <  xt*+l  then 
tna  =  k 
end  if 

continue  k  loop 

Uq  —  N o  Tin  —  p 

loop  from  k  =  X  to  Nv  - 1  {find  index  of  angles  60, 
if  6q  >  Vk  and  90  <  t)*+i  then 
n0  =  k 
end  if 

if  0a  >  Vk  and  0a  <  then 
na  =  k 
end  if 

continue  k  loop 
m  —  ma 
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Fig.  7.  Elevation  dependence  drawing  pattern  for  spherical  coordinates 
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Fig.  8.  Power  pattern  of  uniformly  excited  circular  aperture  viewed  from  0o  =  60° 
and  =  30° 
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Fig.  9.  Power  pattern  of  uniformly  excited  circular  aperture  viewed  from  80  =  60° 
and  <)> o  =  30° 
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loop  from  1=1  to  Nu  —  N,  —  1 
to  =  to  —  (— l)1  (/  —  1) 

if  ( m  <  1)  then  {force  index  to  be  periodic} 
m  =  to  +  Nu  —  1 
else  if  (to  >  Nu  -  1)  then 
to  =  to  —  (Nu  —  1) 
end  if 

if  |ttm  -  <j) 0|  <  7r/2  or  |tim  -  $o|  >  37t/2  then 
i  =  Nv,j  =  0,  L\  —  Nv  -  no,  L2  =  n0-  1 

else 

l  “  Mai  J  —  Tla  lj  L\  —  Tlfl  1  j  L2  ““  iVj  Tig 
end  if 

loop  from  =  1  to  ii 
n  =  i  —  k 

{  get  ( um,vn )  and  other  3  corners  } 

{  convert  to  rectangular  coordinates  -  Table  2  } 
{  scale,  rotate,  and  project  using  Eqn.  (13)} 

{  fill  quadrilateral,  then  draw  perimeter  } 
continue  k  loop 
loop  from  k  =  1  to  L2 
n  =  j  +  k 

{  get  {um,vn)  and  other  3  corners  } 

{  convert  to  rectangular  coordinates  -  Table  2  } 
{  scale,  rotate,  and  project  using  Eqn.  (13)} 

{  fill  quadrilateral,  then  draw  perimeter  } 
continue  k  1c  >p 

if  l  >  Nu  -  N,  -  3  and  JV,  >  0  then 
if  um  >  (j) o  then 

if  (um  >  ‘po  and  um  <  <f)a)  then 
p  —  m 

else 

p  =  TO  +  1 

end  if 

else 

if  (um  >  (j>a  and  um  <  <f> o)  then 
p  =  to  +  1 

else 

p  =  TO 

end  if 

loop  from  q  =  1  to  Nv 
{  form  polygon  with  vertices  ( Up,vq )  } 
continue  q  loop 

{  fill  polygon,  then  draw  perimeter  } 
end  if 

continue  l  loop 
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VIII.  Conclusion 


An  algorithm  has  been  presented  that  allows  the  rapid  plotting  of  an: 
tenna  and  scattering  patterns  in  three  dimensions.  Because  of  the  special 
properties  of  the  types  of  functions  considered,  the  plotting  speed  is  essen¬ 
tially  the  same  as  if  no  hidden  lines  were  removed  Any  graphics  system  that 
allows  a  polygon  fill  operation  may  implement  the  algorithm.  Although  not 
suitable  for  mechanical  pen  plotters,  the  algorithm  is  ideal  for  video  dis¬ 
plays  and  laser  printers.  The  specific  variations  in  rectangular,  cylindrical, 
and  spherical  coordinates  have  been  developed  and  tested  yielding  excellent 
results. 
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Appendix:  Fortran  77  Programs 


This  appendix  contains  the  Fortran  77  programs  RECT3D,  CYLA3D,  CYLR3D, 
and  SPHR3D  as  well  as  all  associated  subroutines.  The  required  inputs  to 
each  program  are  listed  in  the  comments  found  in  the  source  code.  These 
programs  were  used  to  generate  all  the  results  presented  in  this  report.  The 
output  of  each  program  is  a  PostScript  file.  Information  on  programming  in 
the  PostScript  language  is  available  from  most  bookstores. 
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1 

2  C 

3  C 

4  C 

5  C 
e  C 
7  C 
e  C 
9  C 

10  C 

11  C 

12  C 

13  C 

14  C 

16  C 
18  C 

17  C 

18  C 

19  C 

20  C 

21  C 

22  C 

23  C 

24  C 

26  C 
28  C 

27  C 

28  C 

29  C 

30  C 

31  C 

32  C 

33 

34 
36 

36 

37 

38 

39  C 

40  C 

41  C 

42 

43 

44 

46 
40 

47  C 

48  C 

49  C 


PROGRAM  RECT3D 

a********************************************** ******************** 


*  THIS  PROGRAM  PRODUCES  A  POSTSCRIPT  FILE  REPRESENTING  A  THREE  * 

*  DIMENSIONAL  PLOT  OF  A  FUNCTION  IN  RECTANGULAR  COORDINATES.  * 

******************************************************************* 

*  TIMOTHY  J.  PETERS  LAST  UPDATED  * 

*  THE  AEROSPACE  CORPORATION  3/1/91  * 

*  2360  EAST  EL  SEGUNDO  BOULEVARD.  *  * 

*  EL  SEGUNDO,  CA  90246  * 

******************************************************************* 

*  INPUTS:  * 

*  * 

*  NU  -  NUMBER  OF  U  POINTS.  * 

*  NV  -  NUMBER  OF  V  POINTS.  * 

*  AX  -  X  DIRECTION  SCALE  FACTOR.  * 

*  AY  -  Y  DIRECTION  SCALE  FACTOR.  * 

*  AZ  -  Z  DIRECTION  SCALE  FACTOR.  * 

*  PO  -  PHI  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=2PI.  * 

*  TO  -  THETA  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=PI.  * 

*  IC  -  IF  IC=1  THEN  CONVERT  THE  FUNCTION  VALUES  TO  DB.  * 

*  NOTE  THAT  IF  IC=1  THEN  F  MUST  BE  IN  THE  RANGE  0<=F<=1.  * 

*  ETA  -  PLOT  FLOOR  IN  DB.  * 

*  NUU  -  EFFECTIVE  ZERO  IN  DB.  * 

*  U(MAX)  -  U  COORDINATE  ARRAY.  * 

*  V(MAX)  -  V  COORDINATE  ARRAY.  * 

*  F(MAX.MAX)  -  FUNCTION  VALUE  MATRIX.  * 

*  * 

*  OUTPUT:  * 

*  * 

*  RECT.PS  -  POSTSCRIPT  FILE  REPRESENTING  THE  PLOT.  * 

*  * 


******************************************************************* 
„  PARAMETER  (MAX=100) 

REAL+4  U(MAX) ,V(MAX) ,F(MAX,MAI) ,NUU 
OPEN ( UNIT=2 , FILE= ’ RECT . PS ’ ) 

REWIND (2) 

RAD= . 17463293E-01 
PI= . 3141693E+01 

******************************************************************* 

*  READ  THE  INPUTS.  * 

^***t*******************>r****** ******* ***************************** 

OPEN (UNIT=1 , FILE= 'DATA » ) 

READ ( 1 , * )  NU , NV , AX , AY , AZ , PO , TO , IC , ETA , NUU 
READ(1,*)  (U(M),M=i,NU) 

READ(i,*)  (V(N),N=i,NV) 

READ(1,*)  ((F(M,N) ,M=1,NU) ,N=1,NV) 

*  *  *  *  *  *  *  *  *  *  *  *jfc  *  *  ***  *  +  *  *  * + *  4c  *  *  *  %  *  *  *  ***  *  *  *  ***  *  **  *  **  *  *  *  *  ***  *  *  *  *  *  *  **  +  * * ★ 

*  DEFINE  SOME  MACROS  IN  POSTSCRIPT.  * 

******************************************************************* 
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so  WRITE(2,100)  ’initgraplu.es  erasepage  latter’ 

51  WRITE(2,i00)  ’/m  {moveto}  def  /I  {liheto}  del  /s  {stroke}  def’ 

52  WRITE(2,100)  ’/sg  {setgray}  def  /c  {closepath  1  sg  fill  s}  def* 

53  WRITE(2f100)  ’ / w  {closepath  0  sg  s}  def  /t  {translate}  def’ 

64  C  *****************#*******************e*********iic***iM>****"iMi******** 

65  C  *  SET  THE  LINE  CHARACTERISTICS.  * 

no  c  ******************************************************************$ 

67  WRITE(2,i00)  ’0.5  setlinewidth  1  setlinecap  1  setlinejoin’ 

so  C  *  TRANSLATE  THE  ORIGIN  TO  THE  GEOMETRIC  CENTER  OF  THE  PAPER.  * 

60  C  ******************************************************************* 

01  XQFFSET=306 .0 

02  Y0FFSET=396.0 

03  WRITE(2,10i)  XOFFSET.YOFFSET, ’  t' 

64  C  ** %♦*****♦******♦♦*♦***♦****♦♦****+♦***♦****♦♦♦**♦♦*♦** ************ 

05  C  *  CONVERT  THE  OBSERVATION  ANGLES  TO  RADIANS.  * 

06  C  ******************************************************************* 

07  PQBS=RAD*PO 

08  TOBS=RAD*TO 

09  C  ******************************************************************* 

70  C  *  COMPOTE  THE  ROTATION  ANGLES  WHICH  YIELD  THE  DESIRED  OBSERVATION  * 

71  C  *  ANGLES.  * 

72  C  ******************************************************************* 

73  RP=-P0BS-PI/2 . 0 

74  RT=TOBS 

76  CP=COS(RP) 

70  SP=SIN(RP) 

77  CT=COS(RT) 

78  ST=SIN(RT) 

79  C  ******************************************************************* 

so  C  *  IF  REQUESTED  CONVERT  THE  DATA  TO  DB  SCALE.  * 

81  C  *************************  **************  ******  I********************** 

82  IF  (IC  .EQ.  1)  THEN 

83  TR=iO .0**(0 . i*ETA) 

84  DO  i  M=1,N0 

86  DO  2  N=i,NV 

80  IF(F(M,N)  . LE.  TR)  THEN 

87  F(M,N)=(ETA+ABS(NOIJ))/ABS(NOII) 

88  ELSE 

80  F(M,N)=(10.0*AL0G10(F(M,N))+ABS(NUU))/ABS(NUU) 

00  END  IF 

01  2  CONTINUE 

02  1  CONTINUE 

03  ELSE 

04  END  IF 

05  C  ******************************************************************* 

00  C  *  SET  SOME  CONSTANTS.  * 

07  C  ******************************************************************* 

08  IF  ((POBS  .GE.  0.0)  .AND.  (POBS  .LT.  PI/2.0))  THEN 
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89 

100 

101 

102 

103 

104 

105 
100 

107 

108 
1C« 
110 
111 

112  C 

113  C 

114  C 

115 
110 

117 

118 

119 

120 
121 
122 

123 

124 

125 
120 

127 

128 

129 

130 

131 

132 

133 

134 

135 

136 

137 

138  C 

139  C 

140  C 

141 

142 

143 

144  C 

145  C 
148  C 
147 


IS1=1 

IS2=i 

ELSE  IF  ((FOBS  .GE.  PI/2.0)  .AND.  (POBS  .LT.  PI))  THEN 
TSi=-l 
IS2=1 

ELSE  IF  ((POBS  .GE.  PI)  .AND.  (POBS  .LT.  3.0*PI/2.0))  THEN 
ISi=-l 
IS2=-i 

ELSE  IF  ((POBS  .GE.  3.0+PI/2.0)  .AND.  (POBS  .LE.  2.0+PI))  THEN 
ISl=i 
IS2=-1 
ELSE 
END  IF 

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * * *  *  *  *  *  *  <•  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * 
*  BEGIN  SEQUENCE.  * 

******************************************************************* 
DO  3  L=i,NV-i 

IF  ((POBS  .OF,.  0.0)  .AND.  (POBS  .LT.  PI/2.0))  THEN 
MO=L-i 
NO=L 
Mi=L 
Ni=L 

ELSE  IF  ((POBS  .GE.  PI/2.0)  .AND.  (POBS  .LT.  PI))  THEN 
M0=NU-L+1 
NO=L 
M1=NU-L 
N1=L 

ELSE  IF  ((POBS  .GE.  PI)  .AND.  (POBS  .LT.  3.0+PI/2.0))  THEN 
M0=NU-L+1 
N0=NV-L 
M1=NU-L 
N1=NV-L 

ELSE  IF  ((POBS  .GE.  3.0+PI/2.0)  .AND.  (POBS  .LE.  2.0+PI))  THEN 
MO=L-i 
NO=NV-L 
Mi=L 
Ni=NV-L 
ELSE 
END  IF 

***************************************************************** 
*  LOOP  THROUGH  THE  U  VALUES  WITH  V  CONSTANT.  * 

***************************************************************** 
DO  4  K=i,NU-L 
M=M0+IS1*K 
N=NO 

*************************************************************** 
*  COMPUTE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

*************************************************************** 
XS1=U(M) 
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148 

149 
160 
151 
162 
153 
164 
166 

156 

157 
168 
159  C 

100  C 

101  C 

102 

103 
164 
106 
108  C 

107  C 

108  C 
100 

170  C 

171  C 

172  C 

173 

174  4 

176  C 
170  C 

177  C 

178 

179 

180 

181  C 

182  C 

183  C 

184 
,  185 

180 

187 

188 

189 

190 

191 

192 

193 

104 
195 
198  C 


YS1=Y(N) 

ZS1=F(M,N) 

XS2=U(M+1) 

YS2=Vi.N) 

ZS2=F(M+1,N) 

XS3=U(M+i) 

YS3=V(N+1) 

ZS3=F(M+i,N+i) 

XS4=U(M) 

YS4=V(N+i) 

ZS4=F(M,N+1) 

ip*******************^************************  ***************** 

*  ROTATE  THE  4  VERTICES  OF  THE  QUADRILATERAL .  * 

*************************************************************** 

CALL  ROTATK(XSl,YSi,ZSl,Xi,Yi,AX,AT,AZ,RP,RT) 

CALL  R0TAT5(XS2,YS2fZS2,X2,Y2,AX,AY,AZ,RP.,RT) 

CALL  R0TATE(XS3 , YS3 , ZS3 , X3 , Y3 , AI , AT , AZ ,RP ,RT) 

CALL  R0TATE(XS4 , TS4 , ZS4 , X4 , Y4 , AX , AY , AZ ,RP ,RT) 
*************************************************************** 

*  FILL  THE  QUADRILATERAL.  .  + 

*************************************************************** 

WRITE(2,200)  Xl,Yi,X2,Y2fX3,Y3,X4,Y4 
*************************************************************** 

*  DRAW  PERIMETER  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

WRITE(2,201)  X1,T1,X2,Y2,X3,Y3,X4,Y4 
CONTINUE 

***************************************************************** 
*  LOOP  THROUGH  THE  V  VALUES  WITH  U  CONSTANT.  * 

******************************************* *********************** 
DO  E  K=1,NV-L-1 
M=M1 

N=N1+IS2*K 

*************************************  t ************************* 

*  COMPUTE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

**********************************4  **************************** 

XS1=U(M) 

TS1=V(N) 

ZS1=F(M,N) 

XS2=U(M+1) 

YS2=V(N) 

ZS2=F(M+i,N) 

XS3=U(M+1) 

YS3=V(N+l) 

ZS3=F(M+i,N+i) 

XS4=U(M) 

TS4=V(N+i) 

ZS4=F(M,N+i) 

*************************************************************** 
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107 

C 

108 

C 

100 

200 

201 

202 

203 

C 

204 

C 

205 

C 

208 

207 

C 

208 

c 

200 

c 

210 

211 

S 

212 

3 

213 

c 

214 

c 

215 

c 

216 

217 

c 

218 

c 

210 

c 

220 

100 

221 

101 

222 

200 

223 

224 

201 

225 

228 

227 

c 

228 

220 

230 

231 

232 

*  ROTATE  THE  4  VERTICES  OF  THE  QUADRILATERAL .  * 

♦t************************************************* ************ 

CALL  ROTATE(XSl ,TS1 ,ZS1 ,X1 ,Yi ,AX, AT,AZ,RP,RT) 

CALL  R0TATE(XS2,TS2,ZS2,I2,T2,AX,AT,AZ,RP,RT) 

CALL  ROTATE (XS3 ,TS3 ,ZS3 ,X3 ,T3 ,AI,AT,AZ ,RP ,RT) 

CALL  R0TATE(XS4 , TS4 , ZS4 ,14 ,T4 , AX , AT , AZ ,RP ,RT j 
*************************************************************** 

*  FILL  THE  QUADRILATERAL.  * 

*************************************************************** 

WRITE(2,200)  XI ,T1,X2,T2,I3,T3,X4,T4 
*************************************************************** 

*  DRAW  PERIMETER  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

WRITE(2,20i)  X1,T1,X2,T21X3,T3,X4,T4 
CONTINUE 
CONTINUE 

********************************* ********************************** 

*  SHOW  THE  PAGE.  * 

*******************1-1*********************************************** 

WRITE(2,100)  ’ showpage ’ 

******************************************************************* 

*  FORMATS.  * 

******************************************************************* 

F0RMAT(A72) 

F0RMAT(F7 . 2 , II , F7 . 2 , AS8 ) 

F0RMAT(F7.2,1I,F7.2,»  m  > ,F7.2,1X,F7.2, '  1  >,F7.2,1I, 

AF7 .  2 ,  ’  1  ’ ,F7.2,1X,F7.2, ’  1  c») 

F0RMAT(F7 -2,1X,F7 .2, ’  m  * ,F7.2,1X,F7.2, *  1  \F7.2,1I, 

AF7 . 2 , ’  1  \F7.2,1X,F7.2,»  I  w>) 

END 

SUBROUTINE  ROTATE(IA,TA,ZA,X,T,AX,AT,AZ,RP,RT) 

X=AX*COS (RP ) *IA“AT*SIN(RP ) *TA 

T=COS(RT)*(AX*SIN(RP)*IA+AT*COS(RP)*TA)+AZ*SIN(RT)*ZA 

RETURN 

END 
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1  PROGRAM  CYLA3D 

2  C  ******************************************************************* 


3  C  *  THIS  PROGRAM  PRODUCES  A  POSTSCRIPT  FILE  REPRESENTING  A  THREE  * 

4  C  *  DIMENSIONAL  PLOT  OF  A  FUNCTION  IN  CYLINDRICAL  COORDINATES.  * 

5  C  ******************************************************************* 

e  C  *  TIMOTHY  J.  PETERS  LAST  UPDATED  * 

7  C  *  THE  AEROSPACE  CORPORATION  3/1/91  * 

8  C  *  2360  EAST  EL  SEGUNDO  BOULEVARD.  * 

o  C  *  EL  SEGUNDO,  CA  90246  * 

10  C  ******************************************************************* 

11  C  *  INPUTS:  * 

12  C  *  * 

13  C  *  NU  -  NUMBER  OF  U  POINTS.  * 

14  C  *  NV  -  NUMBER  OF  V  POINTS.  * 

16  C  *  AX  -  X  DIRECTION  SCALE  FACTOR.  * 

16  C  *  AY  -  Y  DIRECTION  SCALE  FACTOR.  * 

17  C  *  AZ  -  Z  DIRECTION  SCALE  FACTOR.  * 

18  C  *  PO  -  PHI  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=2PI.  * 

19  C  *  TO  -  THETA  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=PI.  * 

20  C  *  IC  -  IF  IC=1  THEN  CONVERT  THE  FUNCTION  VALUES  TO  DB.  * 

21  C  *  NOTE  THAT  IF  IC=1  THEN  F  MUST  BE  IN  THE  RANGE  0<=F<=1.  * 

22  C  *  ETA  -  PLOT  FLOOR  IN  DB.  * 

23  C  *  NUU  -  EFFECTIVE  ZERO  IK  DB.  * 

24  C  *  NS  -  NUMBER  OF  SEGMENTS  TO  REMOVE.  * 

25  C  *  U(MAX)  -  U  COORDINATE  ARRAY.  * 

20  C  *  V(MAI)  -  V  COORDINATE  ARRAY.  * 

27  C  *  F(MAI.MAX)  -  FUNCTION  VALUE  MATRIX.  * 

28  C  *  * 

29  C  *  OUTPUT:  * 

30  C  *  * 

31  C  *  CYLA.PS  -  POSTSCRIPT  FILE  REPRESENTING  THE  PLOT.  * 

32  C  *  * 


33  C  ******************************************************************* 

34  PARAMETER  (MAX=100) 

36  REAL*4  U(MAX) , V(MAX) ,F(MAX,MAX) ,NUU 

30  0PEN(UNIT=2 ,FILE= ’ CYLA . PS ’ ) 

37  REWIND(2) 

38  RAD= . 17463293E-01 

39  PI=.3141693E+01 

40  C  *♦+******♦**♦♦*♦♦****♦♦♦+*****#**♦#*♦*♦♦♦♦♦**♦**♦♦♦♦♦♦*♦*♦♦**♦♦*♦** 

41  C  *  read  THE  INPUTS.  * 

42  C  ****************************** ************************************* 

43  OPENtUNn^i.FILE^DATA') 

44  READ(1,*)  NU,NV,AX,AY,AZ,PO,TO,IC,ETA,NUU,NS 

45  READ(1,*)  (U(M) ,M=1,NU) 

40  READ(1,*)  (V(N) ,N=1,NV) 

47  READ(1,*)  ((F(M,N),M=1,NU),N=1,NV) 

48  C  ♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ ♦♦♦♦♦♦♦♦♦ 

49  C  *  DEFINE  SOME  MACROS  IN  POSTSCRIPT.  * 
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so  C  I******************************************************************* 
si  WRITE(2,100)  '  initgraphics  erasepage  letter* 

52  WRITE(2,100)  */m  {moveto}  def  /I  {line to}  def  /j  {stroke}  def* 

63  WRITE(2,i00)  ’/sg  {setgray}  del  /c  {closepath  i  sg  -fill  s}  del* 

54  WRITE(2,100)  */w  {olosepath  0  sg  s}  del  /t  {translate}  del* 

55  C  ********************************************-*********************** 

56  C  *. SET  THE  LINE  CHARACTERISTICS.  * 

57  C  ******************************************************************* 

58  WRITE(2,100)  *0.5  setlinesidth  i  set line cap  1  setlinejoin' 

59  C  ******************************************************************* 

80  C  *  TRANSLATE  THE  ORIGIN  TO  THE  GEOMETRIC  CENTER  OF  THE  PAGE.  * 

61  C  ******************************************************************* 

82  XOFFSET=30fl.O 

83  Y0FFSET=390 .0 

64  WRITE(2,10i)  XOFFSET.YOFFSET, *  t» 

85  C  ******************************************************************* 

88  C  *  CONVERT  THE  OBSERVATION  ANGLES  TO  RADIANS.  * 

67  C  ******************************************************************* 

68  PDBS=RAD*PO 

69  T0BS=RAD*T0 

70  C  ******************************************************************* 

71  C  *  COMPUTE  THE  ROTATION  ANGLES  WHICH  YIELD  THE  DESIRED  OBSERVATION  * 

72  C  *  ANGLES.  * 

73  C  ******************************************************************* 

74  RP=-P0BS-PI/2.0 

75  RT=T0BS 

78  CP=COS(RP) 

77  SP=SIN(RP) 

78  CT=COS(RT) 

79  ST=SIN(RT) 

80  C  ******************************************************************* 

81  C  *  IF  REQUESTED  CONVERT  THE  DATA  TO  DB  SCALE.  * 

82  C  ******************************************************************* 

83  IF  (IC  .EQ.  1)  THEN 

84  TR=10.0**(0.i*ETA) 

85  DO  1  M=1,NU 

88  DO  2  N=1,NV 

87  IF(F(M,N)  .LE.  TR)  THEN 

88  F(M,N)=(ETA+ABS(NUU))/ABS(NUU) 

89  ELSE 

90  F(M,N)=(iO.0*AL0G10(F(M,N))+ABS(NUU))/ABS(NUU) 

ot  END  IF 

92  2  CONTINUE 

93  i  CONTINUE 

94  ELSE 

95  END  IF 

98  C  ******************************************************************* 

97  C  *  DETERMINE  THE  VALUE  WHICH  IS  JUST  BELOW  P0+180  DEGREES.  * 

98  C  ******************************************************************* 
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99 

100 

101 

102 

103 

104 

105 
100 

107 

108 

109  3 

110  C 

111  C 

112  C 

113 

114 
116 
110 

117 

118 

119 

120 
121 
122 

123 

124 

125 
128 

127 

128 

129 

130 

131 

132 

133  C 

134  C 

135  C 
138 

137 

138 

139 

140 

141 

142 

143 

144 

145 
140 
147 


IF  (FOBS  .LT.  PI)  THEM 
PA=POBS+PI 
ELSE 

PA=POBS-PI 
END  IF 

DO  3  M=i,NU-l 

IF  ((?A  .GE.  U(K))  .AND.  (PA  .LT.  U(H+i)))  THEN 
IREF=M 
ELSE 
END  IF 
CONTINUE 

******************************************************************* 
*  SEQUENCE  THROUGH  THE  INDICES.  * 

******************************************************************* 

M=IREF 

DO  4  L=1,NU-NS-1 
IS=(-1)**L 
M=M-IS*(L-1) 

IF  (M  .LT.  1)  THEN 
M=M+HU-1 

ELSE  IF  (H  .GT.  NU-i)  THEN 
H=H-(NU-i) 

ELSE 
END  IF 

IF  ((ABS(U(M)-POBS)  .LT.  PI/2.0) 

A  .OR.  (ABS(U(M)-POBS)  .GT.  3.0*PI/2.0))  THEN 

IS=1 
N0=0 
ELSE 
IS=-1 
NO=NV 
END  IF 

DO  B  K=1,NT-1 
N=NO+IStK 

*************************************************************** 
*  GENERATE  THE  RECTANGULAR  POINTS.  * 

*************************************************************** 

cpi=cos(u(h)) 

SPI=5IN(U(M)) 

CPIi=COS(U(M+l)) 

SPI1=SIN(U(H+I) ) 

XSla?(N'i*CPI 

TSi=7(N)*SPI 

ZSi=F(M,N) 

XS2=V(N)*CPI1 

TS2=V(H)*SPI1 

ZS2=F(M+1,N) 

XS3=?(N+1)*CPI1 

TS3=7(N+1)*SPI1 
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148 

140 

160 

161 

162  C 

153  C 

154  C 

155 

156 
1  7 
156 
150  C 

160  C 

161  C 

162 

163  C 

104  C 

105  C 
100 

107  5 

-8  4 

100  C 

170  C 

171  C 

172 

173  C 

174  C 

176  C 
170  100 

177  101 

178  200 
170 

180  201 
181 
182 

183  C 

184 

185 
180 

187 

188 


ZS3=F(M4-1,N4-1) 

XS4=V(N+1)*CPI 

YS4=V(N4-1)*SPI 

ZS4=F(M,N+1) 

****************  I^********************************************** 

*  ROTATE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

CALL  R0TATE(XS1 ,YS1,ZS1 ,X1 ,Yi , AX, AY, AZ ,RP,RT) 

CALL  ROTATE (XS2 , YS2 , ZS2 , X2 , Y2 , AX , AY , AZ , RP , RT) 

CALf  R0TATE(XS3,YS3,ZS3,X3,Y3,AI,AY,AZ,RP,RT) 

CALL  ROTATE (XS4 , YS4 , ZS4 , X4 , Y4 , AX , AY , AZ , RP , RT) 
*************************************************************** 

*  FILL  THE  QUADRILATERAL.  * 

******************************* ************>m* ******** ********* 

WRITE (2,200)  X 1 , Y 1 , X2 , Y2 , X3 , Y3 , X4 , Y4 
*************************************************************** 

*  DRAW  PERIMETER  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

WRITE(2,201)  X1,Y1,X2,Y2)X3,Y3,X4,Y4 
CONTINUE 
CONTINUE 

**************** ****************** *********** ********************** 

*  SHOW  THE  PAGE.  * 

,t  **************************************************************** 

WRITE(2,100)  ’ showpage ’ 

******************************************************************* 

*  FORMATS .  * 

******************************************************************* 

F0RMAT(A72) 

FORMAT (F7 . 2 , IX , F7 . 2 , AB8 ) 

F0RMAT(F7 . 2 , IX , F7 . 2 , '  m  ' ,F7.2,1X,F7.2, '  1  ',F7.2,1X, 

AF7 . 2 , ’  1  1 ,F7.2,1X,F7.2, ’  1  c') 

F0RMAT(F7 ,2,1X,F7 .2,  ’  m  * ,F7.2,1X,F7.2, »  1  ',F7.2,1X, 
ftF7 .2 , ’  1  >,F7.2,1X,F7.2,»  1  w>) 

END 

SUBROUTINE  ROTATE ( XA , YA , Z A , X , Y , AX , AY , AZ , RP , RT) 
X=*X*COS(RP)*XA-AY*SIN(RP)*YA 

Y=C0S(RT)*(AX*SIN(RP)*XA4-AY*C0S(RP)*YA)+AZ*SIN(RT)*ZA 

RETURN 

END 
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1  PROGRAM  CYLR3D 

2  C  ****)|t*){?***  +  ^***)f 

3  C  *  THIS  PROGRAM  PRODUCES  A  POSTSCRIPT  FILE  REPRESENTING  A  THREE  * 

4  C  *  DIMENSIONAL  PLOT  OF  A  FUNCTION  IN  CYLINDRICAL  COORDINATES.  * 

5  c  ******************************************************************* 

0  C  *  TIMOTHY  J.  PETERS  LAST  UPDATED  * 

7  C  *  THE  AEROSPACE  CORPORATION  3/1/91  * 

8  C  *  23S0  EAST  EL  SEGCNDO  BOULEVARD.  * 

9  C  *  EL  SEG'i'UJO ,  CA  9024S  * 

io  C  ******************************************************************* 

n  C  *  INPUTS:  * 

n  C  *  * 

13  C  *  NU  NUMBER  OF  U  POINTS.  * 

14  C  *  NV  -  NUMBER  OF  V  POINTS.  * 

15  C  *  AX  -  X  DIRECTION  SCALE  FACTOR.  * 

io  C  *  AY  -  Y  DIRECTION  SCALE  FACTOR.  * 

17  C  *  AZ  -  Z  DIRECTION  SCALE  FACTOR.  * 

18  C  *  PO  -  PHI  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=2PI.  * 

10  C  *  TO  -  THETA  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=PI.  * 

20  C  *  IC  -  IF  IC=1  THEN  CONVERT  THE  FUNCTION  VALUES  TO  DB.  * 

21  C  *  NOTE  THAT  IF  IC=1  THEN  F  MUST  BE  IN  THE  RANGE  0<=F<=1.  * 

22  C  y  ETA  -  PLOT  FLOOR  IN  DB.  * 

23  C  *  NUU  -  EFFECTIVE  ZERO  IN  DB.  * 

24  C  *  NS  -  NUMBER  OF  SEGMENTS  TO  REMOVE.  * 

25  C  *  U(MAX)  -  U  COORDINATE  ARRAY.  * 

20  C  *  V(MAX)  -  V  COORDINATE  ARRAY.  * 

27  C  *  F( MAX, MAX)  -  FUNCTION  VALUE  MATRIX.  * 

28  C  *  * 

20  C  *  OUTPUT:  * 

30  C  *  * 

31  C  *  CYLR.PS  -  POSTSCRIPT  FILE  REPRESENTING  THE  PLOT.  * 

32  C  *  * 

33  C  ******************************************************************* 

34  PARAMETER  (MAX=100) 

35  REAL*4  U(MAX) , V(MAX) ,F(MAX,MAX) ,NUU 

30  0PEN(UNIT=2,FILE=’ CYLR.PS’) 

37  REWIND (2) 

38  RAD= . 174S3293E-01 

30  PI= . 3141593E+01 

40  C  ************ ******************************************************* 

41  C  *  READ  THE  INPUTS.  * 

42  C  ******************************************************************* 

43  0PEN(UNIT=1,FILE=’DATA’ ) 

44  READ(1,*)  HU , N V , AX , AY , AZ , PO , TO , IC , ETA , NUU , NS 

46  READ(1,*)  (U(M) ,M=1 ,NU) 

40  READ(1,*)  (V(N) ,N=1 ,NV) 

47  READ(1,*)  ((F(M,N) ,M=1 ,NU) ,N=1 ,NV) 

48  C  ******************************************************************* 

40  C  *  DEFINE  SOME  MACROS  IN  POSTSCRIPT.  * 
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BO  C  ******************************************************************* 

61  WRITE(2,100)  ’ initgraphics  erasepage  letter* 

62  WRITE(2, 100)  */m  {moveto}  def  /I  {lineto}  def  /s  {stroke}  def’ 

53  WRITE(2,100)  */sg  {setgray}  def  /c  {closepath  1  sg  fill  s>  def’ 

54  WRITE(2,100)  */w  {closepath  0  sg  s}  def  /t  {translate}  def’ 

55  C  ******************* ******* ******* ft**************** ***************** 

bo  C  *  SET  THE  LINE  CHARACTERISTICS.  * 

57  C  ******************************************************************* 

58  WRITE(2,100)  ’0.6  setlinewidth  1  setlinecap  1  setlinejoin’ 

BO  C  ******************************************************************* 

00  C  *  TRANSLATE  THE  ORIGIN  TO  THE  GEOMETRIC  CENTER  OF  THE  PAGE.  * 

01  C  **.************************+**********+***♦*******+**************** 

02  X0FFSET=308 .0 

63  Y0FFSET=398.0 

04  WRITE(2, 101)  XOFFSET , YOFFSET , ’  t’ 

0B  C  **************************************+**************************** 

00  C  *  CONVERT  THE  OBSERVATION  ANGLES  TO  RADIANS.  * 

07  C  ******************************************************************* 

08  P0BS=RAD*P0 

SO  T0BS=RAD*T0 

70  C  ******************************************************************* 

71  C  *  COMPUTE  THE  ROTATION  ANGLES  WHICH  YIELD  THE  DESIRED  OBSERVATION  * 

72  C  *  ANGLES.  * 

73  C  ***************************************************>•  *************** 

74  RP=-P0BS-PI/2 . 0 

75  RT=TOBS 

70  CP=COS(RP) 

77  SP=SIN(RP) 

78  CT=COS(RT) 

79  ST=SIN(RT) 

80  C  ******************************************************************* 

81  C  *  IF  REQUESTED  CONVERT  THE  DATA  TO  DB  SCALE.  * 

82  C  ******************************************************************* 

83  IF  (IC  .EQ.  1)  THEN 

84  TR=10.0+*(0.1*ETa) 

86  DO  1  M=i ,NU 

86  DO  2  N=1 ,NV 

87  IF(F(M,N)  -LE.  TR)  THEN 

as  F(M,N)=(ETA+ABS(NUU))/ABS(NUU) 

80  ELSE 

oo  F(M,N)=(10.0*AL0G10(F(M,N))+ABS(NUU) )/ABS(NUU) 

91  END  IF 

92  2  CONTINUE 

93  1  CONTINUE 

04  ELSE 

95  END  IF 

90  C  ******************************************************************* 

97  C  *  DETERMINE  THE  VALUE  WHICH  IS  JUST  BELOW  P0+180  DEGREES.  * 

98  C  ******************************************************************* 
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99 

LOO 

101 

102 
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107 

108 

109  3 

110  C 

111  C 

112  C 

113 

114 
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116 
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118 

119 

120 
121 
122 

123 

124 
126  C 

126  C 

127  C 

128 

129 

130 

131 

132 

133 

134 
136 

136 

137 

138 

139 

140 

141 

142 
14  j 
144  C 

146  C 
140  C 

147 


IF  (FOBS  .LT.  PI)  THEN 
PA=POBS+PI 
ELSE 

PA=POBS-PI 
END  IF 

DO  3  M=1,NU-1 

IF  ((PA  .GE.  0(H))  .AND.  (PA  .LT.  0(H+1)))  THF.N 
IREF=M 
ELSE 
END  IF 
CONTINUE 

****************?!************************; ************************** 
*  SEQUENCE  THROUGH  THE  INDICES.  * 

******************************************************************* 
DO  4  K=1 ,NV-1 
N=K 

M=IREF 

DO  S  L=1,NU-NS-1 
IS=(-1)**L 
M=M-IS*(L-1) 

IF  (M  .LT.  1)  THEN 
M=M+NU~i 

ELSE  IF  (H  .GT.  NU-1)  THEN 
H=M-(NU-1) 

ELSE 
END  IF 

*************************%************************************* 

*  GENERATE  THE  RECTANGULAR  POINTS.  * 

*******************  *******>!<**  ***************************  ******* 

CPI=COS(U(M)) 

SPI=SIN(U(M)) 

CPI1=C0S(U(M+1)) 

SPI1=SIN(U(H+1) ) 

XS1=F(H,N)*CPI 

YS1=F(M,N)*SPI 

ZS1=V(N) 

XS2=F(H+i,N)*CPIl 
YS2=F(M+1 ,N)*SPI1 
ZS2=V(N) 

XS3=F(H+1 ,N+1)*CPI1 
YS3=F(M+1 ,N+l)*Sril 
ZS3=V(N+1) 

XS4=F(H,N+1)*CPI 

YS4=F(M,N+1)*SPI 

ZS4=V(N+1) 

***************************************************+*********** 

*  ROTATE  THE  POINTS.  * 

********************************************** ***************** 

CALL  R0TATE(XS1 , YS 1 , ZS1 , XI , Yi , AX , AY , AZ , RP ,RT) 
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148  CALL  ROTATE (XS2 ,  YS2 , ZS2 , X2 , Y2 , AX , AY , AZ , RF , RT) 

146  CALL  R0TATE(XS3 , YS3 , ZS3 , X3 , Y3 , AX, AY , AZ ,RP,RT) 

150  CALL  R0TATE(XS4,YS4,ZS4,X4,Y4,AX,AY,AZ,RP,RT) 

151  C  fc******* t*#****##** ************* ***** **************  t"*********** 

152  C  *  FILL  THE  QUADRILATERAL.  * 

153  C  *<.***********************************************************$* 

164  WRITE(2,200)  XI ,Y1 ,X2 ,Y2 ,X3 , Y3 ,X4, Y4 

165  C  **♦♦**♦***♦***♦*+♦♦♦♦♦+♦♦**♦**♦♦+♦♦*****♦*♦****♦**♦♦****♦♦***** 

1 6tt  C  *  draw  PERIMETER  OF  THE  QUADRILATERAL.  * 

167  C  ******* ♦**♦*****♦♦♦♦♦*****♦******♦ ********************* ******** 

168  WRITE(2,201)  XI ,Y1 ,X2,Y2 ,X3 , Y3 ,X4,Y4 

169  6  CONTINUE 

100  4  CONTINUE 

101  C  ******************************************************************* 

102  C  *  SHOW  THE  PAGE.  * 

103  C  ******************************************************************* 

104  WRITE(2,100)  ’showpage’ 

105  C  ****** ******** ********* ******  I*##* ****** ************* *******  ******** 

100  C  *  FORMATS  * 

107  C  ******************************************************************* 

108  100  F0RMATCA72) 

too  101  FORMAT (F7 . 2 , IX , F7 . 2 , A68 ) 

170  200  F0RMAT(F7.2,1X,F7.2, ’  m  > ,F7.2,1X,F7.2, '  1  »,F7.2,1X, 

171  AF7.2,'  1  \F7.2,1X,F7.2,*  1  c') 

172  201  F0RMAT(F7.2,1X,F7.2, ’  m  ’ ,F7.2,1X,F7.2, '  1  ’.F7.2.1X, 

173  AF7.2, 1  1  \F7.2,1X,F7.2,’  1  u>) 

174  END 

175  C 

170  SUBROUTINE  ROTATE(XA,YA,ZA,X,Y,AX,AY,AZ,RP,RT) 

177  X=AX*COS(RP)*XA-AY*SIN(RP)*YA 

178  Y=COS(RT)*(AX*SIN(RP)*XA+AY*COS(RP)*YA)+AZ*SIN(RT)*ZA 

179  RETURN 

180  END 
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PROGRAM  SPHR3D 

******************************************************************* 


*  THIS  PROGRAM  PRODUCES  A  POSTSCRIPT  FILE  REPRESENTING  A  THREE  * 

*  DIMENSIONAL  PLOT  OF  A  FUNCTION  IN  SPHERICAL  COORDINATES.  * 

******************************************************************* 

*  TIMOTHY  J.  PETERS  LAST  UPDATED  * 

*  THE  AEROSPACE  CORPORATION  3/i/9i  * 

*  2350  EAST  EL  SEGUNDO  BOULEVARD.  * 

*  EL  SEGUNDO,  CA  90245  * 

******************************************************************* 

*  INPUTS :  * 

*  * 

*  NU  -  NUMBER  OF  U  POINTS.  * 

*  NV  -  NUMBER  OF  V  POINTS.  * 

*  U(MAX)  -  U  COORDINATE  ARRAY.  * 

*  V(MAX)  -  V  COORDINATE  ARRAY.  * 

*  F(MAX,MAX)  -  FUNCTION  VALUE  MATRIX.  * 

*  PO  -  PHI  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=2PI.  * 

*  TO  -  THETA  OBSERVATION  ANGLE  IN  RANGE  0<=P0<=PI.  * 

*  IC  -  IF  IC=1  THEN  CONVERT  THE  FUNCTION  VALUES  TO  DB .  * 

*  NOTE  THAT  IF  IC=1  THEN  F  MUST  BE  IN  THE  RANGE  0<=F<=1.  * 

*  ETA  -  PLOT  FLOOR  IN  DB .  * 

*  NUU  -  EFFECTIVE  ZERO  IN  DB.  * 

*  NS  -  NUMBER  OF  SEGMENTS  TO  REMOVE.  * 

♦AX  -  X  DIRECTION  SCALE  FACTOR.  * 

*  AY  -  Y  DIRECTION  SCALE  FACTOR.  * 

*  AZ  -  Z  DIRECTION  SCALE  FACTOR.  * 

*  * 

+  OUTPUT:  * 

*  * 

*  SPHR.PS  -  POSTSCRIPT  FILE  REPRESENTING  THE  PLOT.  * 

*  * 


******************************************************************* 
PARAMETER  (MAX=100) 

REAL*4  U(MAX) , V(MAX) ,F(MAX,MAX) ,NUU 
INTEGER  P.Q 

OPEN (UNIT=2 ,FILE= ’ SPH . PS ' ) 

REWIND (2) 

RAD= . 17453293E-01 
PI= . 3141593E+01 

*  READ  THE  INPUTS.  * 


******************************************************************* 
OPEN (UNIT= 1 ,FILE=’ DATA ’ ) 

READ(1,*)  NU,NV,AX,AY,AZ,PO,TO,IC,ETA,NUU,NS 
READ(1,*)  (U(M),M=1,NU) 

READ(1 ,*)  (V(N),N=1,NV) 

READ(1 ,*)  ((F(M,N) ,M=i,NU) ,N=i ,NV) 

*  *+ *******  *******  *  *  +  sfc  ****  *  ****  *  *****  **  *****************  ******  *****  * 
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50  C  *  DEFINE  SOME  MACROS  IN  POSTSCRIPT.  * 

51  C  ******************************************************************* 

62  WRITE(2,100)  ’initgraphics  erasepage  lotter’ 

53  WRITE(2,100)  ’/m  {moveto}  del  /I  {lineto}  del  /s  {stroke}  del’ 

64  WRITE(2,100)  ’/sg  {setgray}  del  /c  {closepath  1  8g  lill  s}  del’ 

65  WRITE(2,100)  ’/w  {closepath  0  sg  s}  del  /t  {translate}  del’ 

50  C  ******************************************************************* 

67  C  *  SET  THE  LINE  CHARACTERISTICS.  * 

58  C  ******************************************************************* 

59  WRITE(2,100)  ’0.6  setlinewidth  1  setlinecap  1  setlinejoin’ 

00  C  ******************************************************************* 

01  C  *  TRANSLATE  THE  ORIGIN  TO  THE  GEOMETRIC  CENTER  OF  THE  PAPER.  * 

02  C  ******************************************************************* 

03  X0FFSET=308 . 0 

04  YOFFSET=390 .0 

05  WRITE(2,101)  XOFFSET.YOFFSET, ’  t’ 

06  C  ******************************************************************* 

07  C  *  CONVERT  THE  OBSERVATION  ANGLES  TO  RADIANS.  * 

08  C  ******************************************************************* 

09  POBS=RAD*UO 

70  TOBS=RAD*VO 

71  C  ******************************************************************* 

72  C  *  COMPUTE  THE  ROTATION  ANGLES  WHICH  YIELD  THE  DESIRED  OBSERVATION  * 

73  C  *  ANGLES.  * 

74  C  ******************************************************************* 

75  ALPHA=-PI/ 2 . O-POBS 

70  BETA=TOBS 

77  C  ******************************************************************* 

78  C  *  IF  REQUESTED  CONVERT  THE  DATA  TO  DB  SCALE.  * 

79  C  ******************************************************************* 

80  IF  (IC  .Eq.  1)  THEN 

81  TR=10 .0**(0 . 1*ETA) 

82  i/u  I  .."1 1 NU 

83  DO  2  N=1,NV 

84  IF(F(M,N)  .LE.  TR)  THEN 

85  F(M,N)=(ETA+ABS(NUU))/ABS(NUU) 

86  ELSE . 

87  F(M,N)=(10. 0*AL0G10(F(M,N))+ABS(NUU))/ABS(NUU) 

88  END  IF 

89  2  CONTINUE 

oo  1  CONTINUE 

91  ELSE 

92  END  IF 

93  C  ******************************************************************* 

94  C  *  DETERMINE  THE  INDEX  OF  THE  ANGLES  PHI_0  AND  PHI_A  * 

95  C  ******************************************************************* 

90  IF  (POBS  .LT.  PI)  THEN 

97  PA=POBS+PI 

98  ELSE 
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129  C 
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140  C 

141  C 
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143 

144 
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140 
147 


PA=POBS-PI 
END  IF 

DO  3  K=i  ,NtJ-i 

IF  ((POBS  .GE.  U(K) )  .AND.  (POBS  .LB.  U(K+i)))  THEN 
MO=K 
ELSE 
END  IF 

IF  ((PA  .GE.  U(K))  .AND.  (PA  .LT.  U(K+1)))  THEN 
MA=K 
ELSE 
END  IF 
CONTINUE 

*  *  *  *  *  *  *  *  *  *  *  *  **  *  *  **  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  * * * 

*  DETERMINE  THE  INDEX  OF  THE  ANGLE  TOBS.  * 

******************************************************************* 

TA=PI-TOBS 

NO=NV 

NA=HV 

DO  4  K=l,NV-i 

IF  ((TOBS  .GE.  V(K) )  .AND.  (TOBS  .LT.  V(K+1)))  THEN 
NO=K 
ELSE 
END  IF 

IF  ((TA  .GE.  V(K))  .AND.  (TA  .LT.  V(K+i)))  THEN 
NA=K 
ELSE 
END  IF 
CONTINUE 

******************************************************************* 

*  BEGIN  MAIN  LOOP.  * 

****+***+++*+*++++++++++*++** t************************************* 

M=MA 

DO  B  L=1,HU-NS-1 
M=M-((-l)**L)*(L-i) 

IF  (M  .LT.  1)  THEN 
M=M+NU-1 

ELSE  IF  (M  .GT.  HU-1)  THEN 
M=M-(NU-1) 

ELSE 
END  IF 

*  *  *  *  *  *  *  ★  **  *  *  *  *  * + *  *  *  *  *  *  *  **  *  *  *  *  *  *  +  *  *  *  *  ♦  *  *+*  *  *  %  ***  ***  *  +  *  *  +  ★  ***  *  *  ♦  ★  +  * 
*  SET  THE  CONSTANTS  FOR  DRAWING  V.  * 

***************************************************************** 
DP=U(M)-POBS 

IF  ( ( ABS(DP)  .LT.  PI/2.0)  .OR.  (ABS(DP)  .GT.  3.0*PI/2.0))  THEN 
I=NV 
J=0 

L1=NV-N0 

L2=H0-1 
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180 

187 

188 
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100  C 

101  C 

102 

103 
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100  C 


ELSE 
I=NA 
J=NA-1 
L1=NA-1 
L2=NV-NA 
END  IF 

**  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  +  *  *  ****  *  *  *  6  ***#  $  *  *  *  *  *  *  *  *  $  +  *  #***  * 

*  DRAW  THE  V  DEPENDENCY  IN  THE  FIRST  DIRECTION.  * 

♦  ’t'**#*****************************************************  ******* 
DO  6  K=i ,L1 

N=I-K 

*************************** ******************* ***************** 

*  COMPUTE  THE  4  VERTICES  OF  THE  FIRST  QUADRILATERAL.  * 

*************************************************************** 

CPA=COS(D(M)) 

SPA=SIH(a(M)) 

CPB=C0S(U(M+1)) 

SPB=SIN(U(M+1)) 

CQA=COS(V(N)) 

SqA=SIH(V(N)) 

CQB=C0S(V(N+1)) 

SqB=SIN(V(N+l)) 

*  *  *  *  *  *  *  *  *  *  *  **  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  lit  *  *  *  **  *  *  *  *  *  *  *  ]fc  *  *  *  *  *  *  *  *  Hi  *  *  *  *  *  * 

*  COMPUTE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

******** *********** ***** ************************ *************** 

C=F(M,H)*SQA 

XS1=C*CPA 

YSi=C*SPA 

ZS1=F(M,N)*CQA 

C=F(M+l,N)*SqA 

XS2=C*CPB 

YS2=C*SPB 

ZS2=F(M+1 ,H)*CqA 

C=F(M+l,H+l)*SqB 

XS3=C*CPB 

YS3=C*SPB 

ZS3=F(M+i,H+l)*CqB 

C=F(M,H+l)*SqB 

XS4=C*CPA 

YS4=C*SPA 

ZS4=F(M,N+i)*CQB 

*************************************************************** 

*  ROTATE  THE  4  VERTICES  OF  THE  qUADRILATERAL .  * 

*************************************************************** 

CALL  R0TATE(XS1 , YS1 ,ZS1 ,X1 , Y1 , AX, AY , AZ .ALPHA .BETA) 

CALL  ROTATED XS2 , YS2 , ZS2 , X2 , Y2 , AX , AY , AZ , ALPHA , BETA) 

CALL  R0TATE(XS3 , YS3 , ZS3 , X3 , Y3 , AX , AY , AZ , ALPHA .BETA) 

CALL  ROTATE (XS4 , YS4 , ZS4 , X4 , Y4 , AX , AY , AZ , ALPHA , BETA) 
*************************************************************** 
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C 
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C 
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200 

c 

201 

c 
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c 
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c 
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c 
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c 
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c 
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c 
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c 
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c 

222 

c 
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c 
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229 
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232 
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240 

c 

241 

c 

242 

c 

243 

244 

245 

*  FILL  THE  QUADRILATERAL.  * 

**************** ******** ********************************** ***** 

WRITE(2,200)  I1,T1JX2>T2>X3,T3,X4,T4 
*************************************************************** 

*  DRAW  PERIMETER  OF  THE  QUADRILATERAL.  * 

********************************** ***************************** 

WRITE(2 ,201)  X1,Y1,X2,Y2,X3,Y3,X4,Y4 
CONTINUE 

***************************************************************** 
*  DRAW  THE  V  DEPENDENCY  IN  THE  SECOND  DIRECTION.  * 

***************************************************************** 
DO  7  K=1 ,L2 
N=J+K 

*************************************************************** 

*  COMPUTE  THE  4  VERTICES  OF  THE  FIRST  QUADRILATERAL.  * 

*************************************************************** 

CPA=COS(U(M) ) 

SPA=SIN(U(M)) 

CPB=C0S(U(M+1) ) 

SPB=SIN(U(M+1>) 

CQA=COS(V(N)) 

SqA=SIN(V(N)) 

CqB=COS(V(N+l)) 

SQB=SIN(V(N+1)) 

*************************************************************** 

*  COMPUTE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

C=F(M,N)*SQA 

XS1=C*CPA 

YSi=C*SPA 

ZSi=F(M,N)*CQA 

C=F(M+i,N)*SQA 

XS2=C*CPB 

YS2=C*SPB 

ZS2=F(M+1 ,N)*CQA 

C=F(M+1,N+1)*SQB 

XS3=C*CPB 

YS3=C*SPB 

ZS3=F(M+1 ,N+1)*CQB 

C=F(M,N+1)*SQB 

XS4=C*CPA 

YS4=C*SPA 

ZS4=F(M,N+1)*CQB 

*************************************************************** 

*  ROTATE  THE  4  VERTICES  OF  THE  QUADRILATERAL.  * 

*************************************************************** 

CALL  R0TATE(XS1 , YS1 , ZS1 , XI , Y1 , AX , AY , AZ , ALPHA .BETA) 

CALL  R0TATE(XS2 , YS2 , ZS2 , X2 , Y2 , AX , AY , AZ , ALPHA .BETA) 

CALL  R0TATE(XS3 . YS3 , ZS3 , X3 , Y3 , AX , AY , AZ , ALPHA , BETA) 
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CALL  R0TATE(XS4 , YS4 , ZS4 , X4 , 74 , AX , A 7 , AZ , ALPHA , BETA) 

******  **************************  ************  lit****************** 

*  FILL  THE  QUADRILATERAL.  * 

*************************************************************** 

WRITE(2 ,200)  Xi,Yl,X2,Y2,X3,Y3,X4,Y4 

at********************************************************:****** 

+  DRAW  PERIMETER  OF  THE  QUADRILATERAL.  * 

*************************************************************** 
WRITE(2,201)  Xl,Yl,X2tY2,X3,Y3,X4,74 
CONTINUE 

***************************************************************** 
*  CHECK  AND  SEE  IF  A  CUT  IS  REQUESTED.  * 

***************************************************************** 
IF  ((NS  .GT.  0)  .AND.  (L  .GT.  NU-NS-3))  THEN 
*************************************************************** 

*  GENERATE  A  POLYGON  IN  THETA  AT  A  CONSTANT  PHI.  * 

ft******************************  It************************  ******* 

IF  (POBS  .LE.  PI)  THEN 

IF  ( (U(M)  .GE.  POBS) .AND. (U(M)  .LE.  PA))  THEN 
P=M 
ELSE 
P=M+1 
END  IF 
ELSE 

IF  ( (U(M)  .GE.  PA) .AND . (U(M)  .LE.  POBS))  THEN 
P=M+i 
ELSE 
P=M 

END  IF 
END  IF 

CPA=COS(U(P)) 

SPA=SIN(U(P)) 

CQA=COS(V(l)) 

SQA=SIH(V(l)) 

C=F(P,l)*SQA 

XS1=C*CPA 

YS1=C*SPA 

ZS1=F(P,1)*CQA 

CALL  R0TATE(XS1 ,751 ,ZS1 ,X1 , 71, AX, AT, AZ, ALPHA, BETA) 

WRITE(2,300)  11,71 
DO  8  Q=1,NV 
CQA=COS(V(Q) ) 

SQA=SIN(V(Q)) 

C=F(P,Q)*SQA 

XS1=C*CPA 

7S1=C*SPA 

ZS1=F(P,Q)*CQA 

CALL  R0TATE(XS1,7S1,ZS1(XQ,7Q,AX,A7,AZ,ALPHA,BETA) 

WRITE (2, 301)  Xq,7Q 
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332 
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335 
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CONTINUE 

WRITE(2,100)  ’  c  ' 

4**#***>M*** *************** ********************#*********£****% 

*  DRAW  A  LINE  AROUND  THE  PERIMETER  OF  THE  POLYGON.  * 

*************************************************************** 
WRITE(2,100)  *  newpath  > 

WRITE(2,100)  1  0  sg  1 
WRITE(2,300)  XI, Y1 
DO  9  Q=1,NV 
CQA=COS(V(Q)) 
sqA=siN(v(q)) 
c=F(p,q)*sqA 

XS1=C*CPA 

YS1=C*SPA 

zsi=F(p,q)*cqA 

CALL  ROTATE (XSl.YSl, ZS1, Xq.iq, AX, AT, AZ, ALPHA, BETA) 
WRITE(2,30l)  xq.Yq 
CONTINUE 

HRITE(2,100)  ’  s  * 

ELSE 
END  IF 
CONTINUE 

*  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  **  *  *  4C  *  *  *  ****  *  *  *  *  *  *  1)1  *  *  *  *  *  **  *  *  *  **  He  *  *  **  **  *  *  *  *  41  *  *  *  *  4c  *  *  * 

*  SHOW  THE  PAGE.  * 

******************************************************************* 

WRITE(2,100)  ’shonpago1 

******************************************************************* 

*  FORMATS .  * 

******************************************************************* 

F0RMAT(A72) 

F0RMAT(F7 .2 ,1X,F7 .2 , AB8) 

F0RMAT(F7 .2,1X,F7.2, ’  m  ’ ,F7.2,1X,F7.2, * 
fcF7.2, ’  1  ’ ,F7.2,1X,F7.2, ’  1  c») 

F0RMAT(F7.2,1X,F7.2,»  m  * ,F7.2,iX,F7.2, » 
kF7 .2, ’  1  ’ ,F7.2,1X,F7.2, ’  1  *’) 

F0RMAT(F7.2,1X,F7.2, ’  m’) 

F0RMAT(F7.2,1X,F7.2, ’  1») 

END 

SUBROUTINE  ROTATE ( X A , YA , Z A , X , Y , AX , AY , AZ , RP ,RT) 
X=AX*COS(RP)*XA-AT*SIN(RP)*YA 

Y=COS(RT)*(AX*SIH(RP)*XA+AY*COS(RP)*YA)+AZ*SIN(RT)*ZA 
RETURN 
END 
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